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o\ • 

' The instantaneous normal modes corresponding to base pair vibrations (radial modes) and twist 

angle fluctuations (angular modes) of a DNA molecule model at ambient temperature are theoreti- 
cally investigated. Due to thermal disorder, normal modes are not plane waves with a single wave 
number q but have a finite and frequency dependent damping width. The density of modes p(i'), the 
h-^ , average dispersion relation i/{q) as well as the coherence length ^{v) are analytically calculated. The 

■ Gibbs averaged resolvent is computed using a replicated transfer matrix formalism and variational 
' wave functions for the ground and first excited state. Our results for the density of modes are com- 
pared to Raman spectroscopy measurements of the collective modes for DNA in solution and show 
a good agreement with experimental data in the low frequency regime v < 150 cm~^. Radial modes 
extend over frequencies ranging from 50 cm~^ to 110 cm~^. Angular modes, related to helical axis 

^ ■ vibrations are limited to i/ < 25 cm~^. Normal modes are highly disordered and coherent over a 

few base pairs only (^ < 20A) in good agreement with neutron scattering experiments. 
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I. INTRODUCTION 



[ Much attention has been paid in the last twenty years to low frequency dynamics of DNA due to its relevance 
to biological processes. Vibrational modes involving collective motions of groups of atoms have been experimentally 
investigated by means of various techniques as neutron scattering spectroscopy measurements NMR, ... 

Among the latters, Raman studies carried out by H. Urabe et al. [|p0j|] have revealed particulary useful to gain 
^ information on the dependence of low frequency vibrations properties of DNA upon external conditions e.g. water 
content, ionic concentration, temperature. The observation of the Raman scattering intensity in the low frequency 
region is technically very difficult because the strong solvent Rayleigh scattering near zero frequency masks the DNA 
response. Nevertheless a broad band ranging from 60 cm~^ to 100 cm~^ has been evidenced and associated to 
hydrogen bonded base pair vibrations 1^,^. Moreover experiments focusing on oriented solid DNA fibers have also 
exhibit sharp peaks in the Raman intensity at ~ 16 cm"'^, that shifts toward lower frequency region when rising the 
degree of hydration The origin of this peak is far from being obvious. In particular, its interhelical or intrahelical 
' origine has been debated for a long time ||^. 
I To reach a better understanding of the above experimental findings, detailed theoretical analysis of DNA vibrational 
T-H ' motions have been proposed. Devoting a particular attention to hydrogen bond stretching modes, Prohofsky and his 
0^ , collaborators |l^ have been able to find back a vibration mode at 85 cm~^ and confirm its origin. Their approach is 
based on a detailed description of the DNA molecule at the atomic level [0 and a variational calculation method, 
called modified self-consistent phonon approximation (MSPA). Within MSPA, the molecule at ambient temperature is 
reduced to an effective (and complex) harmonic lattice the force constants of which are determined in a self-consistent 
way. This approach has been very useful in calculating many properties of DNA, e.g. the melting of DNA's of different 
sequences , |l^] , the temperature and salt dependence of the B to Z conformation change Q and the stability of 
^ ' drug-DNA complex |jl^. 
Q A qualitative weakness of MSPA, that also arises in other simplified theoretical approaches |]l^ lies in the a priori 

O nature of modes. Though effective elastic constants depends on temperature, normal modes indeed conserve a plane 
^ _ wave structure ]lO| . As a consequence, dispersion relations for the normal modes are well defined as for phonons in 
■ crystals p^ ]. In other words, the coherence length is infinite and the momentum selection rule give rise to a discrete set 
rS I of lines in the theoretical predictions for Raman spectra JTtI . On the contrary, amorphous materials or more generally 
^ ■ thermally disordered systems give rise to normal modes with short coherence lengths, whose power spectrum displays 
a finite width |^,| 18 1. As a result, momentum selection rules break down and light scattering processes may occur 



from essentially all normal modes |l9y2C| ]. Correspondingly for a DNA molecule in solution, the Raman spectrum is 
continuous and does not depend on the wave length of the incident laser beam nor on the scattering angle 

Understanding the vibrations of a system in presence of thermal (or configurational) disorder is a difficult 
task pl| , p2| . Non-harmonicities in the interactions between constituents may modify deeply the usual phonon picture 
of wave propagation in solids and generally forbid any rigorous analytical treatment of the dynamical equations. An 
interesting approach to circumvent this difficulty has been proposed in the context of liquid state dynamics . The 
idea is to start from a randomly chosen configuration at thermodynamical equilibrium. Then, the equations of motion 
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around this fixed initial configuration can be finearized, defining some instantaneous normal modes (INM) and their 
corresponding relaxation times i.e. frequencies. All these quantities may then be averaged over the initial configura- 
tion chosen from the Gibbs ensemble at the desired temperature. One obtains this way a precise description of the 
short times dynamics based on the average spectrum of relaxation times and the statistical properties of the INM . 
As a major advantage, the INM approach takes into account thermal disorder through the choice of the random initial 
configuration and gives rise to highly disordered normal modes with finite autocorrelation length as expected at finite 
temperature. Due to the linearization procedure, the information available through INM calculations is restricted to 
short times dynamics. This limitation is however not serious in the range of frequencies mentioned above for DNA 
collective dynamics as we shall see later. 

From a theoretical point of view, the calculation of the INM around a fixed intial configuration is a hard task that 
can be performed analytically for simple enough models only [p4|-p6[. In this paper we describe the torsional and 
hydrogen bond stretching vibrations of a very simple DNA model by normal modes that are not plane waves using the 
INM approach. We are able in particular to calculate the frequency dependent density of state at a given temperature, 
the coherence length of the normal modes and some pseudo dispersion relations at ambient temperature. The model 
that we consider has been introduced in a previous article to study the DNA denaturation driven by temperature or 
mechanically induced by a torque p?! . This study was motivated by the recent development of micromanipulation 
technique s p^ , |29| that have allowed the direct observation of DNA denaturation induced by an external torsional 
stress p0| , |3l| . Our model describes the DNA molecule at the base pair level, and for each base pair we have just two 
degrees of freedom that is the base pair radius and twist angle. As a result, the potential energy is sufficiently simple 
to allow for sophisticated analytical calculations of normal collective modes. From a technical point of view, our 
calculation is inspired from a recent study of instantaneous modes in an one dimensional disordered system JBS], that 
mixes techniques used in localization theory (the resolvent calculation) , disordered system (replica trick) | |33[ | and a 
variational Gaussian wave function method. As we shall see, the interest of the calculation is two-fold. First, it gives 
theoretical predictions for the spectrum of modes, the effective dispersion relations and damping width at ambient 
temperature that can be directly confronted with Raman spectroscopy and to neutron scattering data. Secondly, it 
is a way to check the validity of (and eventually improve) our model in very different experimental conditions from 
the the ones it was originally designed for. 

The paper is organized as follows. In Section II, we define the DNA model and explain how to compute its 
statistical physics properties. The choice of the force constants and the main thermodynamical features are then 
exposed. Section III is devoted to dynamics and to the definiton of the instantaneous modes. The relationship 
between the latters and Raman spectra is evoked. The analytical framework necessary to calculation the properties 
of INM is presented in Section IV. Results are given and compared to experiments in Section V. Finally, we present 
some conclusions and perspectives in Section VI. 



II. PRESENTATION AND MAIN FEATURES OF THE DNA MODEL 



In this Section, we describe the DNA model under study and briefly recall how it can be studied using statistical 
mechanics techniques as well as its main thermodynamical features, see We then discuss the properties of the 
normal modes, e.g. dispersion relations, spectral density, ... of the molecule in the absence of temperature and 
interaction with the solvent. 



A. Definition of tlie model 



Our model reproduces the Watson-Crick double helix (B-DNA) as schematized fig 1. Each base pair (j = 1, . . . , N) 
is described by its radius rj and the angle (pj in the plane perpendicular to the helical axis . The sugar phosphate 
backbone is made of rigid rods, the distance between adjacent bases on the same strand being fixed to L = 6.95A. 
Conversely, the distance hj between base planes j — 1 and j is allowed to fluctuate. The Hamiltonian V associated to 
a configuration of the degrees of freedom {rj, ipj} is the sum of three different contributions, 

N N 

that we now describe. 



Hydrogen bonds inside a given pair n are taken into account through the short-range Morse potential 35 
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(e-"(^^-^)-l)' (2) 

with R = lOA . The width of the weU amounts to 3a^^ ~ O.SA in agreement with the order of magnitude 
of the relative motion of the hydrogen bonded bases . A base pair with diameter r > Vd = R + 6/a may 
be considered as open. The potential depth D, typically of the order of O.leV jl^ depends on the base pair 
type (Adenine-Thymine (AT) or Guanine-Citosine (GC)) as well as on the ionic strength. Note that the Morse 
potential Vm increases exponentially with decreasing r <. R and may be considered as infinite for r < rmin = 

9.7A 

• The shear force that opposes sliding motion of one base over another in the B-DNA conformation is accounted 
for by the stacking potential 

V;(r,,r,_i) = £;e-^('^^+'--^-2«) (r, - r,_i)' ■ (3) 

Due to the decrease of molecular packing with base pair opening, the shear prefactor is exponentially attenuated 
and becomes negligible beyond a distance ~ 55^^ = 10 A, which coincides with the diameter of a base pair [^ . 

• An elastic energy is introduced to describe the vibrations of the molecule in the B phase, 

V,ir„r,.u0,)^K{h,-H)^ 
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= Kl^^L^- r] - r]_, + 2r,r,.i cos 9, - H) (4) 

where 9j = ipj — ipj^i is the twist angle between base pairs j ^ I and j. The helicoidal structure arises from the 
choice oi H < L: in the rest configuration rj — R at T = OK, Vb is minimum and zero for the twist angle 8 > 
with sin(0/2) — \/L?- — H"^ /{2R). The above definition of Vb holds as long as the argument of the square root 
in (^) is positive, that is if r^, rj_i, 6j are compatible with rigid rods having length L. By imposing Vf, = oo for 
negative arguments, unphysical values of r^, rj_i, 9j are excluded. As the behaviour of a single strand (r > r^) is 
uniquely governed by this rigid rod condition, the model does not only describe vibrations of helicoidal B-DNA 
but is also appropriate for the description of the denaturated phase | |27| |. 

B. Calculation of partition function 

The configurational partition function at inverse temperature (3 = l/lksT) reads 

/■oc /"OO /-oo noo r ^ 

Z= ridri / difi... / riydrN / dipN exp I - 0V[{rj,(pj}][ S{ifi)Y]_x{(^3) ■ (5) 

•^fmin J — OO J V^i^ J — OO L J o — 9 

J — ^ 

The angle of the first extremity of the molecule is set to </5i = with no restriction (due to the arbitrary choice of the 
angular reference axis, see fig 1) whereas the last one ipiq is not constrained. The x factors entering are defined by 
x(6'j) = 1 if < 6*^- = ifj — tpj-i < TT and otherwise to prevent any clockwise twist of the chain. Partition function 
Z can be calculated using the transfer integral method |Q , 

poo 



ndri / TNdrN / dipN {rN,ipN\T |ri,0) , (6) 

n ^ '^min ^ — OO 

where the transfer operator entries read (r, (^|T|r', (^') = T(r, r', 6) with 9 = — ip' and 

r(r, r', 9) = X(r, /) exp {-pVb (r, /, 9)} x{9) , (7) 
X{r,r') = V^ expS^~^{V^{r) + V^{r'))-pV,{r,r')^ . (8) 

At fixed r, r', the angular part of the transfer matrix T is translationally invariant in the angle variables if, f' and 
can be diagonalized through a Fourier transform. Thus, for each Fourier mode k we are left with an effective transfer 
matrix on the radius variables 
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(9) 



with 

Yk{r,r')^ I dOexpi ~ pVt{r,r\9)-ik9\ . (10) 



The only mode contributing to Z is /c = once ipN has been integrated out in (H). In the N ^ oo hmit, the free- 
energy density / is simply given by / = — fesTlnAo, where Aq is the maximal eigenvalue of Tq whose corresponding 
eigenvector will be denoted by '4>o{r). 



C. Thermodynamical properties and parameters 

The ground state wave function ipQ{r) is shown figure 2 at ambient temperature T = 300 K. It is entirely localized 
in the Morse potential well as expected for a DNA molecule in B configuration. The first excited eigenstate of the 
transfer matrix has an excess free-energy AG with respect to ^/^o and is delocalized: it extends over all values of r > r^, 
vanishes for r < and thus represents a denaturated molecule. At some higher temperature T^, the bound wave 
function disappears and ^/jq suddenly undergoes a delocalization transition. In other words, hydrogen bonds break up 
and Tm can be interpreted as the melting temperature 

The values of the parameters entering the potential energy is discussed in [ p7[ and are listed below: 

• inverse hydrogen bond length: a = 6.3A ^ . 

• zero temperature interplane distance: H = 3A. 

• Morse potential depth: D = Q.lGeV. 

• attenuation coefficient for stacking interactions: b = 0.49A ^ . 

• stacking stiffness: E — AeV/A . 

• backbone elasticity constant: K — O.OlAeV/ A . 

The values of a and b have been borrowed from literature |2^. The choice of the other parameters D, E, K, H 
ensures that geometrical and thermodynamical properties as the average twist angle, the mean axial distance between 
successive bases in the B conformation, the melting temperature Tm = 350 K and the denaturation free-energy AG 
are correctly predicted p7| ]. 

Notice that the main uncertainty in this tuning procedure arise with the choice of the stacking stiffness E. Three 



possible pairs of parameters {D, E) that correctly fit Tm — 350 K are listed in Table A 4, as well as the corresponding 
denaturation free-energies AG at T = 300 K. We have selected the pair giving the largest prediction for the denat- 
uration free -ene rgy that is in closest agreement with thermodynamical estimates of AG 139]. It can be easily seen 



from Table A 4 that E varies much more than D and AG and is therefore less accurately predicted than the other 



parameters of the model. 



III. DYNAMICS AND INSTANTANEOUS NORMAL MODES 



Our aim is to perform an analytical calculation of the spectrum of torsional and radial vibrations on characteristic 
time scale of the pico seconds and at ambient temperature. In this section, we first write the dynamical equations for 
the model. We then linearize these equations around a given (and randomly chosen) configuration of the thermally 
equilibrated system and define the instantaneous normal modes as the vibrations of the DNA molecule around this 
configuration. The density of instantaneous normal modes can be related to some extent to the Raman intensity. 
Finally, we consider the special case of zero temperature. The explicit calculation of the dispersion relations and the 
density of modes allows for a decoupling of the angular and radial motions. It will also provides useful comparisons 
with the finite temperature results of Section V. 
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A. Dynamical equations for the DNA model 



We denote the configuration of the molecule at time t by {r*, ip*}. The equations of motion read 



mf* 



1 dV 



where the potential energy V has been defined in ([T]). The effective half mass m = 300 u.m.a. of a base pair 
takes into account the atomic consituents of the nucleotide and of the backbone as well as the primary hydration 
shell which is tightly bound to the base The size of the primary shell depends on the hydration degree 

and is of the order of 10-20 water molecules per nucleotide. The characteristic relaxation time of the primary shell 
is typically ri ~ 10^^*^ s Therefore, for dynamical processes taking place on time scales r < ri, that is for 

frequencies v > 0.3 cm^^ primary shells may be considered as rigidly linked to the bases and simply taken into 
account through the effective mass m. Our choice for m is in close agreement with the mass considered by Volkov 
and Kosevich [Q. Taking into the additional water shell mass and averaging over the possible bases A,T,G and C, 
these authors have calculated somes estimates of the masses of the nucleotide m„ = 199 u.m.a. and of the backbone 
elements m;, ~ 109 u.m.a. giving a total mass m = m„ + mt, = 308 u.m.a. per base |4^ ]. 

The secondary hydration shell contains less rigidly bound water molecules that induce some friction acting on the 
DNA molecule. The characteristics scale time of viscosity dissipation at ambient temperature is T2 ~ 10^^^ s pl| , p3[ , 
that is of the same order of magnitude as in bulk water. The viscous terms —7 f* and —7 (p* that should be included 
in the motion equations (O) can be neglected on time scales t < T2 that is for frequencies v > 30 cm^^ only. In the 



following we shall use the non dissipative equations (11) keeping in mind that the interpretation of the results must 
be made with care for frequencies smaller than 30 cnT^"^. 



B. Linearization approximation and instantaneous modes 



We linearize the equations of motions ( |ll[ ) around a casual initial configuration C = {vj, ipj} of the system already 
in thermodynamical equilibrium, defining for each base pair j, 



= rj + 



(12) 



Once linearized the dynamical equations (|T^) are rewritten in terms of the displacement variables y. 
read 



myj 



R.(j)j and 



(13) 



where 



1 d^V 



2 drjdvk 



1 d^V 



2R drjdipk 



and 



1 d'^V 



V) jk 



2E? dipjdipk 



(14) 



are the N x N matrices of second derivatives of the potential energy V around configuration C. To solve the linear 
system (O) one has to find the eigenvalues A of the 2A^ x 2A^ Hessian matrix 



(15) 



It is important to notice that due to the dependence on the initial configuration C the elements of the matrix T>'^ are not 
translationally invariant. Consequently, the eigenvectors oiT)^ , i.e. the instantaneous normal modes corresponding 
to the initial configuration C are not plane waves. 
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The eigenvalues histogram p^{\) give the density of the normal modes as a function of the initial configuration. 
The latter is distributed according to the Gibbs measure 

V{C = {r„ ^,}) - I exp {~m{rj,Vj}]) , (16) 

where V is the Hamiltonian defined in equation (|^) . Once averaged over distribution (^6|) , the mean density of states 
p(A) is available and depends only on the temperature T = 1//3. The frequency spectrum Pfii') is straightforwardly 
obtained through the relationship, see (|l3|), 

1 



Pf{iy) ^AnV^ p{X) . (17) 

To lighten notation, wc shall drop the / subscript in ( p7| ) and use indifferently the same letter p to denote the density 
of states as a function of the eigenvalue A or frequency ly. 

The spectrum of the Hessian matrix ( p^ ) is not necessarily positive. Some instantaneous modes may be unstable and 
grow with time within the linear approximation. Their corresponding eigenvalues A are negative and the associated 
frequencies v are purely imaginary. Following the conventions of | p3| , imaginary frequencies v will be conveniently 
represented by minus their modulus —{vl, that is by points on the negative frequency semi-axis. We shall come back 
in Section V to these modes. 




C. Relationship with Raman spectra 

In disordered solids the Raman scattering intensity is directly related to the density of normal modes. Indeed when 
the coherence length of the normal modes is short compared to optical wavelengths the conservation of momentum 
is no longer a restrictive selection rule and does not give rise to a discrete set of lines for the spectrum. Light 
scattering processes occur from essentially all the normal modes of the material and the spectra are continuous. As 
we will see later, the assumption of short coherence length is well verified in our model at finite temperature. The 
relationship between the Raman intensity I(V) and the density of normal modes p{i') is generally expressed via the 
light-to- vibrations coupling coefficient Civ) [ [l9| , 

= p{.) C{.) (^^4^) ' 

where 

n{y) = (19) 

is the average population of level v. The spectral dependence of C{v) is still an unsettled question. Early studies 
devoted to the relationship between Raman spectra and densities of modes assumed that the light-to-vibrations 
coupling was indepent of frequency, i.e. C(y) = 1 ||l9[| . Later works conjectured a quadratic behaviour C{v) ~ 
at very low frequencies and a less steep increase for larger u pC| ]. Recently, comparisons between neutron scattering 
experiments and calorimetric measures have give n evidence for a linear dependence C{i') cx v for different glasses in 
the frequency range 8cto~^ < v < 100 cm~^ |44| , |45[| , the upper bound being related to the Debye frequencies of the 
corresponding crystals. 



D. Normal modes at zero temperature 



The linearized equations af motions at zero temperature are obtained performing an expansion of the potential 
energy V up to the second order atound the rest positions: r] = R + y], = n9 + cfi^ R: 



my] = -a^Dy] - Kyy {2y% + y\^^ + y%_^) - E [2y] - y]^^ - y]_^) - Ky^ {4>]^, - cj)]_,) 
m4>] = -K^^ [24)] - (t>]+^ - + Ky^ (y*+i - y]_^) 



(20) 
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where 



Kyy = {KRl/H^) (1 - COS e)' , (21) 
K^^^{KRl/H^){sin^Q^) , (22) 
Ky^ - {KRll H^){8ineo){l - cosQa) (23) 



The plane waves 



exp {i(gn - 27ri/±(g)t)} (24) 



are solutions of (EG) with the relations of dispersion v±{q) showed in fig 3. Due to the difference of the order of 
magnitude between a^D = 6.33 eV A , E = A eV A (or for the other choice of the parameters a?D — 5.93 eV A , 
E = 0.74 eV A~'^) and Ky^ = 18 IQ-^ eV A^^ , Kyy = 6 10"^ eV A'^ , K^^ = 54 10^^ eyA"^ it is clear that the 
angular and radial motions take place on two different time scales and become independent. The dispersion relations 
obtained when setting Ky^ — and Kyy — are indistinguishable from the previous ones and read 



1 a^D ^E ,^ , 

Vr{q) = v+ = —\ + 2-(l-cosg) (25) 

ZTT V TO m 



Ml)^'^--^J^a-cosq) (26) 



The density of states p{u) = 1 / {2tt\i'' {q)\) for each branch is given by 

PrH = , ^""^ (27) 

P^M = ^=^^= ■ (28) 

Note that the above densities are both normalized to 1/2 so that the total density of states Pr{i^) + P^pi^) is properly 
normalized to unity. As shown fig 4, Van Hove singularities are located at frequencies u such that the denominators 



in (28) vanish, that is stationary points of the dispersion relations (M). 



E. Structure of the Hessian matrix 

From the above considerations on the effective decoupling of radial and torsional modes, T)^ will be set to zero 
hereafter. Therefore, the Hessian matrix ( |l5| ) is comprised of two N x N matrices, that is the radial Hessian matrix 
and the angular Hessian matrix P^. In what follows, the notation V'^ will generically refer to either of these two 
matrices. 

Because only nearest neighbors along the molecule interact with each other, the Hessian matrix has a band diagonal 
structure, 

^do{rj,rj^i,0j) + ^do{rj,rj+i,-9j+i) if k = j 



(T)C],.-) -^di{rj,rj^i,9j) if k=j-l , 

^^'■'"^ -H(r,+i,r„0,+i) if + 1 ■ ^^^^ 

if |A:-j|>2 



Explicit expressions of elements c?o ^^nd di are: 

• for the radius Hessian matrix T)^ , the elements 

Mry) = l'^ir) + 2E , 

di{r,r') = 2E , (30) 

do not depend on the twist 9. For simplicity, we have not considered the exponential attenuation term in Vs 
which is almost constant (and equal to unity) in the B-DNA phase. Kyy has been set to zero since it is three 
orders of magnitude smaller than other radial force constants. 
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for the twist Hessian matrix DJ^, 



where the backbone potential has been defined in (^). 



1. ^ 

i?2 de^ 



(r, r',t 



(31) 



IV. ANALYTICAL FRAMEWORK 



A. Definition of spectral quantities 



We call Ag (respectively g) the eigenvalues (respectively the components of the associated eigenvectors normalized 
to unity) of , with e = 1, . . . , iV. Most spectral properties of can be obtained through the calculation of the 
resolvent [BSl 



N C C 



E 



e=l 



A - AC + 



(32) 



where I denotes the identity operator. 
Introducing the density of eigenvalues 



1 ^ 



(33) 



e=l 



we may rewrite the discrete sum over eigenstates e in (|3^) as an integral over eigenvalues with measure . The 
knowledge of the trace of the resolvent, 



1 ^ 1^1 f 

^EGL(A + ..) = ^Ex^IcT7^ = y 

— 1 e— X 



\ — 11 + it ' 



gives then access to the density of states through identity 



p'^(A) = lim Im 



1 ^ 
^EGaa(A + *e) 



Another quantity of interest is the autocorrelation function of eigenvectors at distance d, 

N-d 

C 



(34) 



(35) 



(36) 



We then define A^{\,d) as the average value of Al{d) over all eigenvectors e lying in the range A < A^ < A + dA. 
This autocorrelation function is simply related to the off-diagonal resolvent (|3^) through 



N-d a 
lvEGa,a+.(A + *e)- / 



X — II + ie 



A^{li,d) 



that generalizes equation ( |34| ) to non zero values of d. Taking the imaginary part of equation (pT]), we obtain 

N-d 

p^{\)A^{\d) 



1 



lim Im 

TT £-»0+ 



^ E Ga,a+d(A + «e) 



(37) 



(38) 



Therefore, the calculation of the large distance d behaviour of a+d(A + *e) will give access to the asymptotic scaling 
of the autocorrelation function 



A' 



■(A,d) cx (e""' 



(39) 



and thus to some effective relation of dispersion \'''{q) and to the coherence length 1/(t''(A) of the instantaneous modes 
at finite temperature. 
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B. Average over the instantaneous molecular configuration 



To perform the average over the molecule configurations C, we first rewrite the resolvent ( p2[ ) as the propagator of 
a replicated Gaussian field theory |B3[ 



I Uf=i dxj {xa Xb) exp (A + ie) ^^^^ - f Y^f^^^i '^Ik Xk 



I nf=i dx, exp (i(A + ^e) Ef.i a;/ - § Effe=i 



N 



N 



N 



= l™o"n / Ild^ii^^-^b) exp 2(^ + *^)II^/" 2 ^ ^- 



jk Xj-Xk 



(40) 



j,k=l 



Replicated fields Xj = (xj, . . . , x") are n-dimensional vector fields attached to each site j. The positivity of e ensures 
that the Gaussian integrals in (^) are well defined. 

At equilibrium, the molecular configuration C is drawn according to the Gibbs measure (|l|). We shall denote the 
average of any quantity over distribution (p^ ) in the same way as its configuration dependent counterpart but without 
C subscript. For instance, the average resolvent reads 



GabiX + ie) 



ridri 



dip I 



r^vdr^v / d^N V{{rj,ipj}) G^^'^'''\\ + ^e) 



(41) 



where V is the Gibbs measure (p^. It is important to keep in mind that P^j, vanishes for |j — fc| > 2, see Section 
III. 2. Therefore, only replicated variables Xj,Xk corrsponding to adjacent base pairs along the molecule {j = k± 1) 
interact together in the expression (EO) of the resolvent. 



C. Transfer matrix formalism 



The one dimensional structure of the interactions in (41) can be exploited through the introduction of a transfer 
matrix T relating the molecular variables rj , ipj as well as the replicated variables Xj to their counterparts at site 
j — 1. The entries of this matrix can be read from (|l|,^,^l]). 



{r,^,x\T\r\^\x') ^T{r,r\9) exp [ -{X + ie){x^ + x'^) - -do{r,r\9)x^ 

-'-do{r' ,r,-e)x'^ + ^ di{r,r' ,e)x.x' 



(42) 



In the above expression, do and c?i are respectively the diagonal and off-diagonal elements of V'^ , see ( p9|) while T has 
been defined in (|^). Notice that T is a symmetric (but not Hermitian) matrix. 

By definition, the average resolvent Gat is the correlation function of replicated variables Xa and Xb (^0|). Within 
the transfer matrix formalism, this mean dot product can be computed from the knowledge of eigenvalues A^ and 
all (normalized) eigenvectors '^e{r,(p,x) of T. Calling d the axial distance \a — 6| between sites a and b along the 
molecule, the average resolvent reads in the thermodynamical limit 



GabiX + ie) = lim 



dx dr dip x^ '^'oi^^ 'Pi x) '^e{r, ip, x) 



(43) 



The maximal eigenvalue (in modulus) of T coincides with £ = Q and increasingly excited states correspond to £ > 1. 
The average diagonal resolvent Gaa and thereby the mean density of states may be obtained from (|4^) when d — 0, 



— lim Re 

TT £^0+,n-*0 



J dx dr dip 4'o(?', ip, x)^ {x^)'^ 



(44) 



Expression (^) may also be used to compute the autocorrelation function of eigenvectors at large distances d. for 
which the resolvent (|4^) scale asymptotically as, see (|39|), 



Ga,a+d{X + ie) cx exp(-(T(A) d + i g(A) d) 



(45) 
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where both a{> 0) and q depend on the energy level A. When d is large, the sum in ( p^ is dominated by the £ = 1 
contribution and equation ( [45[ ) may be reformulated as 

lim i In Ga,a+d{>^ + *e) = -cr(A) + i q{X) (46) 

d — *oo (J 

which allows to derive a and q from the knowledge of Aq and Ai. We shall explicitely compute the inverse length a 
and wave number q in Section V and compare them to the zero temperature results of Section II. D. 

Notice that a and q defined in ( p5| ) do not exactly coincide with the values of a'^ and q'^ appearing in (|3^) averaged 
over C with distribution ([l^). To obtain the latter, one should average the logarithm of the resolvent (j' (quenched 
average) and not simply compute the logarithm of the mean resolvent as in (|4^) (annealed average) |Q. Due to 
concavity of the logarithm function, a defined in (^|) is not only a approximation but also a lower bound to the 
average value of . Comparison with numerical simulations made for the so-called "smallworld" lattice have shown 
that a provides a very good estimate of the mean a'^ p2| . To sum up, we can compute from the average resolvent 
G an upper bound l/cr to the coherence length of the eigenvectors at level A as well as an estimate q of their typical 
wave number. 



D. Rayleigh-Ritz formula and variational approach 



The diagonalization of the transfer matrix and the analytical continuation in n — > could in principle be performed 
exactly due to the rotational invariance of T in the n-dimensional space of replicated variables. To avoid this tedious 
calculation, we resort to a Gaussian variational approach whose accuracy and reliability has been recently validated 
in the case of diffusion on random lattices p5| , p2| . 

We start from the Rayleigh-Ritz formula for the largest eigenvalue Aq of a (real-valued) matrix T, 



An = max7?.(*) 



where 



7^(4') 



A lower bound Ag to Aq can be obtained through maximization of the r.h.s. of 
wave functions 5'g((3) parametrized by some tunable variables Q. According to 



m 

(48; 



(48) 



(49) 

within a suitable trial family of 
E9I), we obtain 



where the optimal parameter Q 



opt 



Ag=7^(vI/g(Q„p,)) 
is the solution of 



dQ 



(50) 



(51) 



Notice that within formulation (|5C|,|5l|), T needs no longer to be real- valued and can assume complex values. 
Hereafter, we apply this procedure to the transfer matrix (p2|) and resort to the following Ansatz, 



tpoir) exp [-Qir) 



(52) 



where V'o('') is the maximal wave function of the Tg transfer matrix (^) of Section II. B. This choice ensures that the 
correct eigenvalue Aq is recovered when the number of replicas n strictly equals zero: Ao(n) = Aq + Oin). The x 
dependence of the trial wave function (^2|) results from the similarity of the operator T (^ ) with the transfer matrix 
of a ID chain of interacting spherical spins, for which the ground state wave function is Gaussian pO{ . Furthermore, 
this Ansatz allows for an exact calculation of the n-dimensional integral over x and an easy continuation of the result 
to real values of n. Note here that the wave function ( p2| ) does not depend on Lp. This is no assumption for n = 0, 
see Section II. B and supposed to give quantitatively good results for n > 0. 
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The calculation of the Rayleigh-Ritz functional TZ ( |49| ) with the Gaussian Ansatz (^2|) is exposed in Appendix A. 
The resuhing functional optimization equation over Q reads 

- ' dr de T[r,r ,0) -tpolr ) — — , (53) 



2Q(r) Jo ' ' ' W(r,r',0) 

where 

W(r, r', d) = (g(r) + A + ^e - do(r, r', 0))(g(r') + A + ze - do(7-', r, ~Q)) - di(r, r', 0)^ . (54) 
The solution of the above equation gives access to the variational estimate of the density of modes from (p^, 



p'W^--J ^'^^oW ■ (55) 

We can pursue the procedure to compute the excited eigenvalues A^ of T. Drawing our inspiration from the ID 
spherical spins chain transfer matrix |40|| , we look for a variational wave function of the type 

*,^(r,^,x)-*g(r,(^,f)P/(x\x2,...,x") , (56) 



where is a polynomial of n variables of total degree £. To obtain the asymptotic behaviour of the resolvent (|4E|), 
we only need the first excited state ^ = 1, whose corresponding variational polynomial is clearly Pf{x) = x^, giving 
thereby from ( |49| ) 

r d. (r r' 

^-} = -2 I / drdr' / d9 T{r,r',9) Mr) Mr') ' ' 



K Jo ' ' ' ^"'^^ ^""^ ^ W(r,r',0). 

x|ao^" drVo(0'QM"'} (57) 
in the n — > limit as shown in Appendix A. 

V. RESULTS AND COMPARISON WITH EXPERIMENTS 
A. Numerical resolution of the self-consistent equations 

We first compute the ground state wave function ipoi'r) and the corresponding eigenvalue Aq, see Section II. B, by 
means of Kellog's method The integration range [rminirmax] over r is discretized into a set of points r^, 

a = 1, . . . , with ri = rmin = 9.8A and r„^ = rmax = 10.7X. 

Self-consistent equations ( |53| ) for Q{r) can then be solved iteratively. As can be checked at zero temperature, 
convergence is improved by iterating the equations for 1/Q{r) rather than for Q{r) itself. We stop the iteration 
process as soon as the differences between the I /Q{r a)^s and their images through the iteration become smaller than 
lO"'^ for aU 1 < a < n^. 

Numerical difficulties come from the limits e — > and — > cxd. This can be best seen for base pairs vibrations in 
the simplest case E — 0. The exact solution to equations ( |5^ ) is Q(r) = A + ie — ^V^{r). In other words, eigenvalues 
A and radii r are in simple correspondance : to any permitted eigenvalue A is associated one (or a few) radius r(A) 
such that A — ^Vj/^{r{X)) and the density of states reads 

p{\) = lim - dr t^J'^^ = HTT^TtW Mr{X))^ ■ (58) 



In practice however, the integral in (Bq) is discretized as follows 



j_ Y> Mrg)^ e 



p(A) = iimzy: .. rr.;.. • (59) 



Consider the eigenvalue A^ corresponding to radius for some arbitrary integer 7. Dominant 0(1/ e) contributions 
to the density yo(A^) in (B3) come from a in the range 7— A<a<7 + A with A ~ 2 e/V"'{ry). Problems arise 
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when A is close to unity. When A scans the interval [A-y; A-y+i], the index of the only (for A = 1) contributing term 
to p(A) jumps from a = 7toa = 7 + lat some intermediate A which will be a local minimum of the density p. Such 
local fluctuations are pure artifacts of the discretization procedure and must be removed by keeping A ^ 1, that is 
l/rir <C e <C 1. Typical suitable values of the parameters are e = 2.10~^, — 4.10'^, giving a spectrum p{\) almost 
normalized to unity (with a small error ~ e). 

Note that the same reasoning holds for the numerical calculation of the torsional spectrum. The discretization of 
the integral over Q < 9 < -k must be replaced by a sum involving a large number ng = 900 of terms to reach a good 
accuracy of the results. 

B. Radial modes 

Radial mode spectra obtained for different choices of D and E are displayed figure 5. They exhibit smooth shapes 
and Van Hove divergences have been smeared out by thermal disorder, compare to figure 4. As at zero temperature, 
the overall width of the spectrum is an increasing function of E. 

For E — 4eV/ A , the general form of p{h') is reminiscent of the density of states at zero temperature, with a shoulder 
in ~ 75 cm^^ and a maximum in ~ 135 cm^^, in correspondence with the edges of the zero temperature 
spectrum, i^*^ — 73.8 cm^^ and u'^ = 138.7 cm^^. A careful analysis even show a quantitative agreement between the 
density of states at T = 300 K and T = K for frequencies lying in the range 110 cm~^ < ly < 130 cm~^. A similar 
behaviour, that is the robustness of the central part of the spectrum to (weak) disorder was also observed in [ p2[ . 

o 2 — 

At a weaker stiffness E = O.lAeV/A , a single bump is observed. The range of allowed frequencies at zero temper- 
atures ly'^ = 71.5 cm~^ < 1/ < 87.5 cm~^ is indeed too narrow and both peaks merge under the action of thermal 
disorder. Note that a very small fraction of modes seem to be unstable and give rise to imaginary frequencies, see 
Section III.B. We however discard them since their integrated sum is smaller than the accuracy e of the calculation. 

Figure 6 shows the dispersion relations at ambient temperature for radial modes. The frequency i^riQ) is an 
increasing function of the wave number over the interval < q < tt. Since the range of allowed frequencies is much 
larger at T = 300 K than at T = K, there is no general coincidence with the zero temperature dispersion relations as 

„ 2 

can be seen for E — OJAeV/ A . For larger stiffness constants, the dispersion relations for both temperatures however 
coincide for medium wave numbers, i.e. when 1 < q < 2 roughly. In agreement with the above analysis of the density 
of states, the effective thermal disorder gets weaker and weaker as the stiffness constant E grows. 

When E increases, the molecule becomes more and more rigid since radii Vj less and less fluctuate from base pair 
to base pair. The wave function ipo{r) gets more and more concentrated around the minimum of the Morse potential 
r = R, see figure 2 and the region of integration over r that mostly contributes to the density of states in ( |55| ) becomes 
narrower and narrower. On the opposite, for small E, ipo^r) mainly reflects the structure of the Morse well whose 
flanks are not accessible at zero temperature. The tails of ■00 are large and give rise to some tails for the density of 
states. The cross-over between both regimes takes place dX E D a? , that is of the order of a few eVj A . 

The inverse coherence length cr^ is plotted figure 7 as a function of frequency. In the central ragion of the spectra, 
the corresponding autocorrelation lengths are ^ ~ 0.7 for E ~ 0.74 and ^ ~ 4 for E — A with a more sensitive 
dependence on v in the latter case. As expected from the above discussion, f increases with E. The values of the 
coherence length ^ are in good agreement with the equilibrium correlation distance Q defined through 

((r,-(r))(r,+,-(r)))(xe-''/« (d ^ oo) . (60) 

A thermodynamical calculation of Q can be carried out from the knowledge of the excited states of the transfer matrix 
To confined to the Morse potential Results are ( 0.51 for E = 0.74 and C 0.93 for = 4 Note that C 
is the inverse damping with of the static structure factor whereas ^(i') is an energy (frequency) dependent coherence 
length. 

C. Angular modes 

We now turn to the angular spectrum. The density of modes is shown on figure 8. We first concentrate on positive, 
that is real frequencies. The band edge = 13.6 cm~^ of the zero temperature spectrum visible on figure 4 disappears 
at finite temperature. A sharp maximum now takes place at vm — 7 cm^^. As expected from zero temperature 
calculations, the width of the peak is smaller than for radial modes and can be estimated to /S.v ~ 15 cm~^. 

The densities of states a.t E — 0.74 and E — A coincide within 0.1%. We have numerically checked that the angular 
mode spectrum depends extremely weakly on the stiffness constant E over the whole range Q < E < oo. In other 
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words, unlike radial modes, angular modes are not sensitive to the width of the ground state wave function i/jq, i.e. 
to fluctuations of the base pair radius r, see figure 2. This observation is supported by inspection of the variational 
parameter Q{r) entering wave function (^2|). At fixed frequency, both imaginary and real parts of Q{r) are indeed 
almost constant on the whole range of radius 9.8A < r < 10.7A. 

The robustness of the spectrum displayed fig. 8 can be understood by looking at the variations of dalr, r' , 9) around 
the thermal average positions r = r' — (rn) and 9 = {9n). The second derivative of do with respect to r at fixed 

twist angle equals d^do — 0.0'i5eV/A while the range of fluctuations of r is given by the largest squared width of 

ipo (corresponding to i? = 0) and reads A(r^) ~ 0.04A . Therefore the fluctuations of the radius r modifies do by 

^9^do.A(r^) ~ O.OOleV/ A , that is by less than 2% of do typically. Repeating the same calculation for twist-induced 

fluctuations, we obtain dgdo ~ 20eV/A and A(6'fj) ~ ksT / [d^R^) ~ 0.005 rad^ . The resulting variations of do due 

to changes of twist are of the order of 5i9^do.A(6'^) ~ Q.QheV/ that is comparable to do. 

Consequently, an excellent approximation of the angular mode spectrum may be obtained by the following simple 
argument. Let us call 

M6)^^'^{R.R,o) (61) 

the diagonal element of T)'^ (|2^) where we have for simplicity identified (r„) with R since (r-„) — i? ~ O.OlA ^ i? at 
T = 300 K). The twist angle 9 is approximately distributed with the Gibbs measure, see (|lO|), 



p{0) = Y„{R,R) ^ ' ^' ^) r ■ (62) 



Then we may substitute the variational equation (p3) on Q{r) with 



r M Q + X + ie-do{9) 



which involves a single parameter Q. Solving equation (|63|), the density of angular modes equals —lm{l/Q)/n and is 
in excellent agreement with figure 8. From a numerical point of view, this approximation is much less time consuming 
than the full resolution of (^3|). In fact, equation (^) required the integral over 9 to be computed for each value of r 
(and at each iteration step) and the solving time was therefore rir times larger than for (|63|). 

Close to zero frequency, the density of states vanishes as Pi^(i^) c>c since the density of eigenvalues p(A) is finite 
in A = 0, see (p7|). We now turn to negative, that is imaginary frequencies. Computing the integrated density of 
unstable modes, we see that the latters represent roughly 20% of angular modes. They extend down to frequencies 
equal to v- = —20 cm"-'^ with a maximum in ~ —3.5 cm^^. We shall come back to the physical implications of 
these modes in next Section. 

The relation of dispersion Vip{q) for the angular modes is displayed figure 9 over the whole range of real and 
imaginary frequencies. We also show on figure 10 the inverse autocorrelation length aip{v). As for radial modes, 
the dispersion relations at T = 300 K gets close to its zero temperature counterpart at intermediate wave numbers 
0.5 < q < 1.5 corresponding to frequencies cm~^ < < 10 cm"^. This coincidence is accompagnied by small values 
of a on this interval of frequencies, giving rise to a coherence length of the order of ^ ~ 2. Conversely, large frequencies 
correspond to highly disordered modes: q > 2 and monotonously increasing cr, with ^ ~ 0.4 at = 30 cm~^. Notice 
that the statical correlation length is the same as for radial modes, see Section V.B. 

Unstable modes have also short autocorrelation lengths, e.g. ^ ~ 0.4 at i/ = —20 cm^^. However, their wave 
numbers are much smaller and can be considered as contant (and zero) for < —10. The autocorrelation function of 



unstable eigenmodes (36) therefore do not change sign over a typical distance dg ~ Tr/q ^ 1. Unstable modes can be 
seen as unstable acoustic phonons, involving coherent rotations of the molecule extending over regions of size dg. 

To end with, notice that the dispersion relation and the inverse coherence length both exhibit an inflection point 
in 1/ = as shown fig. 9 and 10. This is an artefact of the representation of unstable modes as negative frequencies. 
Consider for instance the damping width a close to zero frequency. In the natural A eigenvalue parametrization, we 
expect a non singular behaviour in the vicinity of A = 0: cr(A) = ctq + <7iX + O(A^). For positive eigenvalues A, the 
frequency u is defined as oc VA (|i^ while using the negative-imaginary convention of Section III.B, i/ oc —\/—X for 
negative eigenvalues. The expansion of the inverse coherence length as a function of (small) frequencies thus reads: 
ai^iy) = (To + d'iiy.\h'\ + 0(|t/p) and is singular with an inflection point in i/ = 0. 
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D. Discussion 



Our calculation shows that the coherence length of normal modes is finite and remains of the order of unity. 
Eigenvectors at finite temperature are thus far from being plane waves as in the zero temperature case. The power 
spectrum, that is the Fourier transform of the autocorrelation function at frequency v ( p5| ) acquires a Lorentzian 
form centered around a certain wave number q(y) with a width oiy). This behaviour is experimentally observed 
in neutron scattering experiments |p],|l8|. Furthermore, the calculation justifies a 'posteriori the absence of selection 
rule on momentum in Raman experiments. Indeed, the coherence length of the disordered phonons is of the order of 
< 20 A and is neghgible with respect to optical wavelengths ~ 4800 A, see Section III.C. 

Both radial and angular spectra are superposed in figure 11. The total spectrum p^v) equals Pri^^) + Pip{v) is 
normalized to unity (half of modes originate from angular degree of freedom, and the remaining half from radial 
ones). The torsion peak appears much more narrow and higher than its radial counterpart. The width of the latter 

amounts to Av ~ 60 cm~^ (respectively lS.v ~ 100 cm~^) for E — 0.74eV/A (resp. for E — AeV/A ) whereas 
(stable) torsional modes spread over a range of Ai/ ~ 15 cm~^. While for E = AeV/ A , torsional and radial spectra 
do not intersect each other, there is an overlap region at smaller stiffness constant around ~ 30 cm~^. Nevertheless, 
both fluctuations take place on basically two distinct scale times. Angular vibrations are associated to a typical 
frequency of ~ 10 cm~^, or equivalently to a typical time ~ i 10^^^ s. Radial modes are present at frequencies 
~ 100 cm~^ that is involve dynamical processes on a scale ~ 0.3 lO"^'^ s. At low frequencies \v\ < 25 cm~^, the 
presence of unstable modes threatens the validity of the linearization procedure used in Section III.B. A more refined 
analysis taking into account non linear terms in the dynamical equations is needed to understand how unstable modes 
(with low wave numbers q) are coupled to stable mode (with larger q) and modify the frequencies of the latters. the 
influence of viscous friction upon modes in this frequency regime would also deserve to be further studied as mentioned 
Section III. A. Conversely, the absence of unstable modes at frequencies larger than i' > 25 cm~^ indicate that the INM 
predictions are reliable for time scale smaller than 10~^^ s. This is enough to identify the range of allowed frequencies 
for torsional modes: v < 30 cm~^. The INM prediction for the location of the torsional peak z/m — 7 cm~^ cannot 
be trust blindly but it reasonably lies in the middle of the zero temperature band, see Section III.D. 

We now turn to the comparison with spectroscopy measurements. To establish the link between the density of states 
and the Raman intensity, the knowledge of the light-to-vibrations coupling C{i') is necessary as seen Section III.C. In 
the absence of precise information on the latter, we have rescaled the density of modes yo(z/) according to formula ( p^ ) 
for the three different choices C(i') = 1, C(i') = ly and C{i') — (note that the second hypothesis for C is the most 
plausible one). The resulting theoretical Raman intensities are shown figure 12. a for E = 0.74eV/A and figure 12. b 

for E = AeV/ A . From the one hand, the overall shape of the spectra change with C with respect to p. In particular 
for C{i') — 1 and C(i^) = i^, the torsional peak diverges at small frequencies and the radial modes bump acquires a 
shoulder from, located in the right flank of the torsional modes. On the other hand, the band of allowed frequencies 
remains unaltered by the choice of C and extends over 50 < < 110 cm^^ for E = 0.74 and 50 < < 150 cm~^ for 
i? = 4. In the latter case, the radial region of the spectrum exhibit two local maxima at the same height for C(z^) = v. 
The rescaled spectrum of fig 12. a for C(i^) = v closely agrees with intensity curves obtained from Raman experiments 
at T = 50°C and shown figure 2 of Ref. |^ on the range of frequencies 50 < < 110 cm^^. The latter measurement 
was performed on calf thymus in solution (10 mM PHB, ph=7). This allows us to think that E = 0.74 is a better 
choice for the stiffness constant than E = A . 

Our results do not show drastic variations with temperature for T ranging from 300 K to the melting temperature 
Tfn = 350 K. At Tm the double helix structure disappears and so do all angular and radial vibrations. Since 
denaturation is described as a first order phase transition by the present model, there is no gradual destabilization of 
the modes and no relevant change in p can be seen as mentioned above. The experimental Raman measurements shown 
on figure 2 of Ref. nevertheless show a smooth change in the shape of the intensity curves over the temperature 
interval 50°C< T < 80°C. This apparent paradox can be easily explained by the fact that our model describes a 
homogeneous sequence of bases. For the latter, the fraction of open base pairs versus temperature exhibits an abrupt 
jump from zero to one at |3^. In the case of a heterogeneous sequence as in the experiments of Ref. ||], the 
denaturation temperature of AT rich regions can be inferred to be Tm = 50°C \ 39 1 and the fraction of open base pairs 
smoothly increases from zero to one over the range of temperatures 50°C< T < 77°C. 

We now turn to angular modes. Raman measurements on fibers at 100% of relative humidity (r.h.) have given 
evidence for a narrow band at 16cm^^, see fig. 1 of Ref. 0|. This peak shifts down to lower frequencies (lOcm^^) 
when increasing the hydration degree of DNA in gels and disappears in the central component for DNA in solution . 
Other measurements on B-DNA fibers al lower ~ 80% r.h. by Lindsay et al. have reported a similar band at a higher 
frequency v ~ 25cm~^ that shifts to lower values as the r.h. increases in good agreement with Urabe et al.'s data [||. 
The observed softening of the frequence may come from the weakening of external e.g. interhelical interactions as 
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well as from the influence of the (rigid) primary and (viscous) secondary water shells, see Section III. A. Note that 
this modes seems to appear when the scattering vector AQ gets parallel to the helical axis. Other molecular system 
as crystals of ATP or guanosine monosphosphate (GMP) that present a columnar stacking of bases (even without 
hydrogen bond interactions) exhibit similar low frequency spectra. The intensity of the low frequency (< 20cm~^) 
mode was found to depend mostly on the stacking degree of bases in a column Our calculation predicts, despite 
the presence of unstable INM and the absence of friction terms in the dynamical equations (11) that angular modes 
are predominant for v < 25cm~^. At such low frequencies, radial modes can be neglected and twist fluctuations are 
responsible for collective vibrations of interplane distances In agreement with experiments, our calculation thus 
predicts that low frequency modes {v < 25cm~^) are related to vibrations of the base pairs column. 

A quantitative comparison of our dispersion relations for angular modes with experimental results is difficult to 
the lack of available data. Neutron scattering measurements have shown that pseudo-dispersion relations with a 
finite damping width can be obtained for low frequencies but these experiments have been performed on crystalline 
DNA fibers. Our model (and the values of parameters exposed Section II. C) are valid for DNA in solution where 
the interaction with surrounding water is radically different and inter helical effects are absent. More precisely, we 
find that at fixed momentum q, the experimental frequency v for fibers is larger than the theoretical predictions for 
diluted DNA. The additional mass due to primary water shells for DNA in solution may account for the reduction 
of frequency to a large extent Conversely, the extrapolated value of a is constant and equal to 0.48 in good 
agreement with the theoretical minimum cr = 0.5 shown fig. 10. 



VI. SUMMARY AND CONCLUSION 



In this article, we have shown how to apply the INM framework to a simple model of DNA molecule and reproduce 
the main features of collective vibration modes at finite temperature. Scale time separation between atomic vibrations 
and collective modes allows to obtain good results at low frequencies without resorting to a detailed description of 
DNA at the atomic level. This simplified modelization of DNA permits in turn a deeper analytical understanding of 
the structure of the modes than e.g. within MSPA. 

The model we had introduced to reproduce thermally and mechanically-induced DNA denaturation transitions has 
proven to be also capable of describing accurately the pico-second dynamics seen through spectroscopy measurements. 
This robustness has been obtained without any modification of the model or any new fit of the constant force or 
geometrical parameters. Remarkably, the comparison with Raman experiments has permitted us to decide the value 
of the only parameter known with some uncertainty from the denaturation experiments, that is the stacking stiffness 
E. The success of the present model to account for completely different experiments comes from its mesoscopic nature, 
that lies at an intermediate level between microscopic modelizations, e.g. studies by Prohofsky et al. flOfl or numerical 
simulations by Lavery and collaborators and elasticity theories as the Worm-Like-Chain model [[49[ and its recent 
extensions pO[ . 

Our calculation gives access to the density of modes p{y) and some statistical properties of the normal modes as the 
dispersion relation v{q) and the autocorrelation length l/a{v). The dispersion relations provide the power spectrum 
of the modes which exhibit a finite damping width a at ambient temperature. 

Let us summarize briefly our main quantitative result. Through a rescaling of the density of modes taking into 
account the light-to-vibrations coupling C, we have related the density of modes to Raman intensity measurements. 

o 2 

The choice of parameters D = QA^eV, E — Q.lAeV/A offers a good agreement with the experiments by Urabe et 
al. (see figure 12. a and figure 2 in ||^). The range of frequencies 50 cm~^ < v < 100 cm~^ corresponding to radial, 
i.e. base pair stretching modes do not qualitatively depend on C. Furthermore, it remains roughly unchanged in the 
whole interval of temperatures 0°C< T < Tm where Tm = 77°C is the denaturation temperature. Indeed, our model 
describes a homogeneous sequence for which the melting transition is very abrupt and not smooth as for disordered 
DNA. 

It would be very interesting to compare our theoretical results for (Jr{v) and Vr{q) with neutron scattering exper- 
iments which to our best knowledge are not available for DNA in solution over the range of frequencies mentioned 
above. A fundamental feature of the modes is that the coherence length ^ is of the order of unity: decorrelation 
between components of the same mode takes place on a few Angstroms. This prediction agrees well with results for 
the static correlation distance C, obtained from statistical mechanics calculations as seen Section V.B. 

As for angular modes, the predicted characteristic frequencies v < 25 cm"-'^ coincide with the Raman measurements 
ly < 16 cm^^ The value of the coherence length in the center of the spectrum, ^ = 1/cr ~ 6.8A is compatible 
with data obtained through neutron scattering experiments on DNA fibers Q. However, as far as angular modes are 
concerned, the present approach suffers from two weaknesses . First, we have not taken into account in the dynamical 
equations the viscous forces that might become relevant at very low frequencies. Secondly, INM can become unstable 
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at small v and the linearization approximation we have used throughout the study reveals dangerous. Further 
information about the non linear couplings between modes would be extremely useful to circumvent this difficulty. It 
is however not clear how such a study could be technically pursued. 
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APPENDIX A: VARIATIONAL CALCULATION AND SELF-CONSISTENT EQUATION FOR Q 



In this Appendix, we compute the Rayleigh-Ritz functional for both the ground state (|5^) and the first excited ( |56| ) 
(with £ — 1) wave functions. We then derive the saddle-points equations on the variance Q{r) for both instantaneous 
modes families. 

1. Calculation of the ground state 

Inserting the Gaussian Ansatz ( |5^ ) in (^9|), we obtain 

drdr' j;;deT{r,r\e) Mr) Mr') {W(r, r', 0)}""/' 

7?. *n = """ 



J^°°, drMr)^ {~8mQ{r)} 



-n/2 



= Ao + |r^[0] + 0(n2) (Al) 

where W is the determinant of the two by two matrix A4 defined by 

f Qir) + X + i€-do{r,r',d) diir,r',9) \ 

di{r,r',e) Q{r') + X + te-do{r',r,-0) J ' 

and 



r'[Q] = Xo I drVoW In Q(r) 

drdr' [ d9 T(r,r' ,0) Mr) Mr') ^■aW{r,r' ,9) (A3) 



up to an irrelevant additional constant. The vanishing condition ( |5l| ) on the functional derivative of r^[(5] with respect 
to Q{r) leads to equation (^3|). 



2. Calculation of the first excited state 

We now compute Af using Ansatz (^6|) and expression (^). The denominator of 7?.[\E'f] reads 

dr J dip J dxMr)^ exp y-^Qir) x 

= iC drMr)^Q{r)-^ , (A4) 

when 71 — > 0. C — dip could be made finite by limiting the range of the twist angle. Such a regularization is 
however not necessary since C also appears in the numerator of ([49), 
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poo pOC 

(^ff|T|*^) = / drdr' dipdifi' T{r,r',Lp-~ if') iPoir) Mr') >^ 

J Vniin J —OO 

X J dxdx' xix[ exp ^ f ' ^ ( f' ) 

drdr' d9 T{r,r' ,0) Mr) M^') (m-'^^ {r,r',e) 



as the number of replicas n vanishes. Using Ag — s- Ao when n ^ 0, we obtain equation (| 



(A5) 



3. Case of base pairs vibrations 



For base pair vibrations, the second derivatives do and di are given in (30) and do not depend on the relative twist 
9 between successive base pairs. Inserting ( pO| ) into the extremization equation (js^) for Q{r) and using the definition 
of the effective radial transfer matrix Tqi we obtain 



Aq Mr) 
2Q(r) 



,,rj., ,,,,,, Q{r')+\ + ie~\V:i,{r')~2E 
dr Mr, r)Mr) ^^^^ 



where 



Wr{r, r') = {Q{r) +\ + ie- ^V:^{r) - 2E^ (^Q(r') + A + *e - ^V^ir') ~2E]~ AE^ 
The ratio (|57|) of the first two eigenvalues of T is given by 



A? 

K 



= -AE 



drdr' ^^^^ Mr) Mr') \ / \ j dr Mr)^Q{r)-' 



(A6) 



(A7) 



(A8) 



4. Case of twist angle fluctuations 

For twist fluctuations, the elements do,di of the matrix of second derivatives depend on radii r,r' as well as 
on the angle 9. No simplification arises as in the base pairs case. The saddle-point equation for Q{r) is precisely 
equation (^) but care must be paid to the backbone potential Vb{r,r' ,9). As can be seen from definition (||), the 
angular integral in (^3) is restricted to twist angles 9 fulfilling the condition 



cost' > 



2rr' 



(A9) 



Indeed, when the inequality in ( [A9| ) is not satisfied, the potential Vf, is infinite, see discussion of Section II. A. The 
same prescription holds for the angular integral in (pT]). 
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TABLES 



D 


E 


AG 


ieV) 


ieV/A^) 


(fcsT) 


0.15 


0.74 


0.762 


0.16 


4 


0.825 


0.17 


12 


0.691 



Table 1: Three choices of the depths of the Morse potential D and of the stacking stiffnesses E giving the desired 
melting temperature = 350 K. The corresponding denaturation free-energies AG are expressed in unit of ksT at 
T = 300 K. 



19 



FIGURE CAPTIONS 



Figure 1: The helicoidal DNA model: each base pair is modehzed through its radius rj and angle ipj. The axial 
distance hj between successive base pairs planes varies while the backbone length along the strands is fixed to L. 

Figure 2: Ground state wave functions ipo{r) at T = 27°C for two choices of the energetic parameters: D — 0.15eT^, 
E = Q.74eV/A^ (full curve) and D = O.WeV, E = AeV/J^ (dotted curve). 

Figure 3: Dispersion relations at zero temperature. From bottom to top: ^^^{q) for torsion modes, Vr{q) for radial 
modes with two different choices of the parameters D = OAbeV, E = O.lAeV/A (full curve) and D = 0.16eV^, 
E = AeV/A (dotted curve). 

Figure 4: Density of states at zero temperature. From left to right: torsion modes spectrum p^{v), radial modes 
spectra p^(i/) for D = 0.15eV, E = Q.7AeV/A^ (fuU curve) and D = O.WeV, E = AeV/A^ (dotted curve). 

Figure 5: Spectra Pr{i^) for radial modes at temperature T — 27°C. Parameters are: D = O.lbeV, E = O.lAeV/ A 

2 

(full curve) and D — QAQeV, E — AeV/A (dotted curve). Each spectrum is normalized to one half. Note the small 
oscillations on the tails of the curves due to numerical artifacts, see Section V.A. 

Figure 6: Dispersion relations i^r{q) for radial modes at temperature T = 27°C. Parameters are: D = 0A5eV, 
E = O.lAeV/A (full curve) and D = 0.16ey, E — AeV/A (dotted curve). For comparison, dispersion relations at 

2 ° 2 

zero temperature are recalled {E = O.TieV/A : dashed curve, E = AeV/A : dash-dotted curve). 

Figure 7: Inverse coherence length cr,.(j/) for radial modes at temperature T = 27°C. Parameters are: D = 0.15eT^, 
E = Q.74eV/A^ (fuU curve) and D = OAGeV, E = AeV/A^ (dotted curve). 

Figure 8: Spectra Pif,{iy) for angular modes for D = 0.16ey, E = ieV/A . The curve for D = 0A5eV, E = O.TieV/A 
is indistinguishable from the latter and from the approximation (^,^2|,^3|). The spectrum is normalized to one half. 
Negative frequencies represent unstable modes according to the convention exposed in Section III.B. 

Figure 9: Dispersion relation i'ip{q) for angular modes (full curve). For comparison, the dispersion relation at zero 
temperature is recalled (dashed curve). Negative frequencies represent unstable modes according to the convention 
exposed in Section III.B. Note the inflection point &t v — due to this representation, see Section V.C. 

Figure 10: Inverse coherence length (J^{v) for angular modes. Negative frequencies represent unstable modes 
according to the convention exposed in Section III.B. Note the inflection point ai v = Q due to this representation, 
see Section V.C. 

o 2 

Figure 11: Spectrum p(i^) for both radial and angular modes Parameters are: D — QAbeV, E — O.TieV/A (full 

o 2 

curve) and D = 0A6eV, E = AeVA (dotted curve). Each spectrum is normalized to unity 

„ 2 

Figure 12: Rescaled spectra for diff'erent light-to-vibrations coupling C for parameters D = 0.15eT^, E = O.lAeV/ A 

a 2 

(a) and D = O.lGeF, E = ieV/A (b). Original spectra p^v) are recalled (full lines) and rescaling functions are: 
C(z^) = 1 (dotted curve), C(z/) — v (dashed curve) and C(y) = (dashed-dotted curve). Vertical units are arbitrary, 
all curves have been multiplied by a constant to meet (and equal unity) in = 50 cm"^. 
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